# ------------------------------------------------------------------------------#
# Description: Create 
# Figure A.6: Heterogeneous peer effects varying the number of households in peer group by
# technology and peer technology
# Figure A.7: Heterogeneous peer effects varying the timing of adoption in peer group by technology
# and peer technology
# Project:     Peer Effects in Voluntary Environmental Policies:
#              An Application to Urban Water Quality
# Author:      Daniel A. Brent | dab320@psu.edu
# Created:     2025-06-11
# ------------------------------------------------------------------------------#

# Clear environment ------------------------------------------------------------
rm(list = ls())
gc()

# Load required packages -------------------------------------------------------
library(pacman)
p_load(data.table, fixest, ggplot2)

options(scipen = 999)

# Load and merge base data -----------------------------------------------------
load("data/adoptions/adoptions.RData")
load("data/bucket/peer_adopt_buckets_any.RData")
dat <- merge(adopt, peer.adopt.bucket.any, by = c("pin", "year"))
rm(adopt, peer.adopt.bucket.any)

load("data/bucket/peer_adopt_buckets_rg.RData")
dat <- merge(dat, peer.adopt.bucket.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.rg)

load("data/bucket/peer_adopt_buckets_any_rg.RData")
dat <- merge(dat, peer.adopt.bucket.any.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg)

load("data/bucket/peer_adopt_buckets_cs.RData")
dat <- merge(dat, peer.adopt.bucket.cs, by = c("pin", "year"))
rm(peer.adopt.bucket.cs)

load("data/bucket/peer_adopt_buckets_both.RData")
dat <- merge(dat, peer.adopt.bucket.both, by = c("pin", "year"))
rm(peer.adopt.bucket.both)

load("data/bucket/peer_elig_buckets_cs.RData")
dat <- merge(dat, peer.buckets.cs, by = c("pin", "year"))
rm(peer.buckets.cs)

load("data/bucket/peer_elig_buckets_rg.RData")
dat <- merge(dat, peer.buckets.rg, by = c("pin", "year"))
rm(peer.buckets.rg)

load("data/bucket/total_peers.RData")
dat <- merge(dat, total.peers, by = "pin")
rm(total.peers)

load("data/parcel_build/parcel_build.RData")
dat <- merge(dat, parcel.build[, .(pin, zip5, sub.area, sqft, lot, beds, baths,
                                   yearbuilt, val.land, val.improve,
                                   ren.pre90, ren.90.99, ren.00.09, ren.10,
                                   water.prob)], by = "pin")
rm(parcel.build)

load("data/census/all_parcels_census_2010.RData")
parcels.2010[, blk := as.numeric(paste0(block, blkgrp))]
dat <- merge(dat, parcels.2010[, .(pin, tract, blkgrp, block, blk)])
rm(parcels.2010)

# Scale peer adoption variables ------------------------------------------------
for (i in c(100, 80, 60, 40, 20)) {
  dat[, paste0('peer.adopt.any.', i, '.scale') := get(paste0('peer.adopt.any.', i)) / 100]
  dat[, paste0('peer.adopt.any.rg.', i, '.scale') := get(paste0('peer.adopt.any.rg.', i)) / 100]
  dat[, paste0('peer.adopt.rg.', i, '.scale') := get(paste0('peer.adopt.rg.', i)) / 100]
  dat[, paste0('peer.adopt.cs.', i, '.scale') := get(paste0('peer.adopt.cs.', i)) / 100]
  dat[, paste0('peer.adopt.both.', i, '.scale') := get(paste0('peer.adopt.both.', i)) / 100]
  dat[, paste0('elig.peers.avg.', i) := get(paste0('elig.peers.cs.', i)) + get(paste0('elig.peers.rg.', i))]
}

dat <- dat[year > 2010]
dat[, elig.peers.avg.100 := elig.peers.cs.100 + elig.peers.rg.100]

# Figure A.6: Heterogeneous peer effects by peer group size ----------------------------------------

for (i in c(100, 80, 60, 40, 20)) {
  for (g in c('any', 'any.rg', 'cs', 'both')) {
    assign(paste0('m.', g, '.', i),
           felm(get(paste0('adopt.', g)) ~ get(paste0('peers.', i)) | blkgrp + year |
                  (get(paste0('peer.adopt.any.', i, '.scale')) ~ get(paste0('elig.peers.avg.', i))) |
                  blkgrp,
                data = dat[elig.cs == 1 & post.rainwise == 0]))
  }
}

load('data/iv/iv_diag_numpeers.RData')

for (g in c('any', 'any.rg', 'cs', 'both')) {
  assign(paste0('coef.plot.', g),
         data.table(
           num = c(20, 40, 60, 80, 100),
           adopt = rep(g, 5),
           beta = sapply(c(20, 40, 60, 80, 100), function(i) get(paste0('m.', g, '.', i))$coefficients[2]),
           ci_low = sapply(c(20, 40, 60, 80, 100), function(i) iv_diag_numpeers[[paste0('m.', g, '.', i)]]$AR$ci[1]),
           ci_hi = sapply(c(20, 40, 60, 80, 100), function(i) iv_diag_numpeers[[paste0('m.', g, '.', i)]]$AR$ci[2])
         ))
}

coef.plot <- rbind(coef.plot.any, coef.plot.any.rg, coef.plot.cs, coef.plot.both)
coef.plot[, Adopt := factor(adopt, levels = c('any', 'any.rg', 'cs', 'both'),
                            labels = c('Any', 'Rain Garden', 'Only Cistern', 'Both'))]

theme_set(theme_bw(base_size = 30))
theme_update(legend.position = 'bottom')

gg <- ggplot(coef.plot, aes(x = num, y = beta)) + 
  geom_point(position = position_dodge(width = 5), size = 3.5, color = 'dodgerblue') +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_hi), linewidth = 2, width = 6.5,
                position = position_dodge(width = 5), color = 'dodgerblue') + 
  facet_grid(~Adopt) +
  xlab("Size of Peer Group") + ylab("Effect on Adoption (percentage points)") +
  scale_x_continuous(limits = c(16, 104), breaks = seq(20, 100, by = 20))

print(gg)

ggsave('output/figures/figure_a6.png', gg, width = 180, height = 150,
       units = "mm", dpi = 150, scale = 2)

# Figure A.7: Heterogeneous peer effects by timing of adoption --------------------------------------

load('data/time/peer_adopt_t_any.RData')
dat <- merge(dat, peer.adopt.t.any, by = c('pin', 'year'))
rm(peer.adopt.t.any)

load('data/time/peer_adopt_t_any_rg.RData')
dat <- merge(dat, peer.adopt.t.any.rg, by = c('pin', 'year'))
rm(peer.adopt.t.any.rg)

load('data/time/peer_adopt_t_cs.RData')
dat <- merge(dat, peer.adopt.t.cs, by = c('pin', 'year'))
rm(peer.adopt.t.cs)

load('data/time/peer_adopt_t_both.RData')
dat <- merge(dat, peer.adopt.t.both, by = c('pin', 'year'))
rm(peer.adopt.t.both)

for (i in 1:5) {
  dat[, paste0('peer.adopt.any.', i, '.scale') := get(paste0('peer.adopt.any.', i)) / 100]
  assign(paste0('m', i),
         felm(adopt.any ~ peers.100 | blkgrp + year |
                (get(paste0('peer.adopt.any.', i, '.scale')) ~ elig.peers.avg.100) |
                blkgrp,
              data = dat[elig.cs == 1 & post.rainwise == 0]))
}

load('data/iv/iv_diag_timing.RData')

for (g in c('any', 'any.rg', 'cs', 'both')) {
  for (i in 1:5) {
    assign(paste0('m.', g, '.', i),
           felm(get(paste0('adopt.', g)) ~ peers.100 | blkgrp + year |
                  (get(paste0('peer.adopt.any.', i, '.scale')) ~ elig.peers.avg.100) |
                  blkgrp,
                data = dat[elig.cs == 1 & post.rainwise == 0]))
  }
}

for (g in c('any', 'any.rg', 'cs', 'both')) {
  assign(paste0('coef.plot.', g),
         data.table(
           num = 1:5,
           adopt = rep(g, 5),
           beta = sapply(1:5, function(i) get(paste0('m.', g, '.', i))$coefficients[2]),
           ci_low = sapply(1:5, function(i) iv_diag_timing[[paste0('m.', g, '.', i)]]$AR$ci[1]),
           ci_hi = sapply(1:5, function(i) iv_diag_timing[[paste0('m.', g, '.', i)]]$AR$ci[2])
         ))
}

coef.plot <- rbind(coef.plot.any, coef.plot.any.rg, coef.plot.cs, coef.plot.both)
coef.plot[, Adopt := factor(adopt, levels = c('any', 'any.rg', 'cs', 'both'),
                            labels = c('Any', 'Rain Garden', 'Only Cistern', 'Both'))]

theme_set(theme_bw(base_size = 30))
theme_update(legend.position = 'bottom')

gg <- ggplot(coef.plot, aes(x = num, y = beta)) + 
  geom_point(position = position_dodge(width = 0.25), size = 5, color = 'dodgerblue') +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_hi), linewidth = 2, width = 0.5,
                position = position_dodge(width = 0.25), color = 'dodgerblue') + 
  facet_grid(~Adopt) +
  xlab("Years since peer adoption") + ylab("Effect on Adoption (percentage points)")

print(gg)

ggsave('output/figures/figure_a7.png', gg, width = 150, height = 150,
       units = "mm", dpi = 150, scale = 2)
